# coding:utf-8
# Author : hiicy redldw ~~ hi ldw
# Date : 2019/02/12
import numpy as np
import scipy
from sklearn import metrics
from sklearn.manifold import MDS
from sklearn.datasets import load_digits


# 计算矩阵各行之间的欧式距离；x矩阵的第i行与y矩阵的第0-j行继续欧式距离计算，构成新矩阵第i行[i0、i1...ij]
def calculate_distance_matrix(x, y):
    d = metrics.pairwise_distances(x, y)
    return d

def main():
    # 距离矩阵D m*m
    data = [[2,1,4,3,5],
         [3,2,6,3,0.8],
         [0.5,3,5,1,6],
         [2,5,0.4,1,5],
         [0.9,2,3,9,2]]
    D = calculate_distance_matrix(data,data)
    D = np.array(D)

    mds = MDS(n_components=3) # FAQ:为什么不一样
    result = mds.fit_transform(D)
    print(f'sklearn MDS算法：',result)
    dpie = 3
    m = numSamples = len(D)
    numFeatures = len(D[0])

    Rlist = [] # 存每一行的
    for i in range(numSamples):
        tdisti = (np.sum(np.square(D[i])))/m
        # print(f'第{i}行的和:',np.sum(np.square(D[i])))
        Rlist.append(tdisti)

    Clist = [] # 存每一列的
    for j in range(numFeatures):
        tdistj = (np.sum(np.square(D[:,j])))/m
        # print(f'第{j}列的和',np.sum(np.square(D[:,j])))
        Clist.append(tdistj)
    td = np.square(D)
    Ddist = np.sum(td)/(m**2)

    DS = np.square(D)

    # 内积矩阵 B
    B = np.zeros((m,m))
    print(len(Rlist),len(Clist))
    for i in range(m):
        for j in range(numFeatures):
            B[i,j] = (DS[i,j] - Rlist[i] - Clist[j] + Ddist)/(-2)
    eigVals, eigVects = np.linalg.eigh(B)
    # eigVectT = eigVects.T
    feaValArr = np.zeros((dpie,dpie))
    assert dpie < m,"降维失败"

    feaindex = np.argsort(eigVals)
    feaindex = feaindex[::-1][:dpie]  # 从大到小
    feaVal = eigVals[feaindex]  # 新特征值 特征向量矩阵
    feaVectarr = eigVects[:,feaindex]  # 特征向量[样本数，特征向量数]
    print(f'feaVectarr:',feaVectarr)
    for i in range(dpie):
        feaValArr[i,i] = feaVal[i]
    feaValArr = np.sqrt(feaValArr)
    Z = np.dot(feaVectarr,feaValArr)
    print(f'Z:',Z,Z.shape)
def other(n=3):
    def calculate_distance(x, y):
        d = np.sqrt(np.sum((x - y) ** 2))
        return d

    def cal_B(D):
        (n1, n2) = D.shape
        DD = np.square(D)  # 矩阵D 所有元素平方
        Di = np.sum(DD, axis=1) / n1  # 计算dist(i.)^2
        Dj = np.sum(DD, axis=0) / n1  # 计算dist(.j)^2

        Dij = np.sum(DD) / (n1 ** 2)  # 计算dist(ij)^2
        B = np.zeros((n1, n1))
        for i in range(n1):
            for j in range(n2):

                B[i, j] = (Dij + DD[i, j] - Di[i] - Dj[j]) / (-2)  # 计算b(ij)
        return B
    data = [[2, 1, 4, 3, 5],
         [3, 2, 6, 3, 0.8],
         [0.5, 3, 5, 1, 6],
         [2, 5, 0.4, 1, 5],
         [0.9, 2, 3, 9, 2]]
    D = calculate_distance_matrix(data,data)
    D = np.array(D)
    B = cal_B(D)
    Be, Bv = np.linalg.eigh(B)  # Be矩阵B的特征值，Bv归一化的特征向量
    print(f'Bv:',Be)
    # print np.sum(B-np.dot(np.dot(Bv,np.diag(Be)),Bv.T))
    Be_sort = np.argsort(-Be)
    Be = Be[Be_sort]  # 特征值从大到小排序
    Bv = Bv[:, Be_sort]  # 归一化特征向量
    Bez = np.diag(Be[0:n])  # 前n个特征值对角矩阵
    # print Bez
    Bvz = Bv[:, 0:n]  # 前n个归一化特征向量
    print("Bvz",Bvz)
    Z = np.dot(np.sqrt(Bez), Bvz.T).T
    print(Z)


if __name__ == "__main__":
    main()
    other()